﻿using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Data;
using System.Drawing;
using System.Linq;
using System.Text;
using System.Windows.Forms;
using DotSpatial.Controls;
using DotSpatial.Topology;
using DotSpatial.Symbology;
using DotSpatial.Data;
using System.Collections;

namespace LiDARAspect
{
    public partial class LiDARAspect : Form
    {
        public DotSpatial.Data.Raster DEMRaster = new DotSpatial.Data.Raster();
        public DotSpatial.Topology.KDTree.KdTree _tree = new DotSpatial.Topology.KDTree.KdTree(2);
        public DotSpatial.Data.FeatureSet LiDARFeatureSet = new DotSpatial.Data.FeatureSet(FeatureType.Point);//LiDAR points or Synthetic Points
        public DotSpatial.Data.Raster GridFile = new DotSpatial.Data.Raster();
        public DotSpatial.Data.FeatureSet OutLet = new DotSpatial.Data.FeatureSet();
        public FeatureSet SyntheticPoints = new FeatureSet(FeatureType.Point);
        public FeatureSet SlopePoint = new FeatureSet(FeatureType.Point);
        public FeatureSet NetLine = new FeatureSet(FeatureType.Line);
        public int numofshape;
        public int nNodes = 800;
        public Node[] Nodes;
        public ArrayList randonList = new ArrayList();
        public LiDARAspect()
        {
            InitializeComponent();
        }

        private void loadDEMToolStripMenuItem_Click(object sender, EventArgs e)
        {
            IMapRasterLayer RLayer = map1.AddRasterLayer();
            GridFile.Open(RLayer.Filename);
            //RLayer.Symbolizer.ShadedRelief.IsUsed = true;
            //RLayer.Symbolizer.ShadedRelief.ElevationFactor = 1;
        }

        private void addOutletToolStripMenuItem_Click(object sender, EventArgs e)
        {
            OpenFileDialog ofd = new OpenFileDialog();
            ofd.Filter = "SHP|*.shp";
            if (ofd.ShowDialog() != DialogResult.OK) return;
            OutLet.Open(ofd.FileName);
            map1.AddLayer(ofd.FileName);            
        }

        private void creatSyntheticPointsToolStripMenuItem_Click(object sender, EventArgs e)
        {
            Random r = new Random();
            ArrayList list = new ArrayList(nNodes);

            for (int i = 0; i <nNodes; i++)
            {
                list.Add(i);
            }

            for (int i = 0; i < nNodes; i++)
            {
                int index = r.Next(0, list.Count);
                randonList.Add(list[index]);
                list.RemoveAt(index);
            }
            int Row = (int)GridFile.Maximum;
            int Col = (int)GridFile.Minimum;

           // DotSpatial.Symbology.RandomExt rnd ;//= new RandomExt();
            nNodes += OutLet.Features.Count;
            
            Nodes = new Node[nNodes];          

            DotSpatial.Data.RcIndex RI = new DotSpatial.Data.RcIndex();

            double[] dbl = new double[2];
            for (int i = 0; i < nNodes; i++)
            {
                Nodes[i] = new Node();
                if (i < OutLet.Features.Count)
                {
                    Nodes[i].X = OutLet.Features[i].Coordinates[0].X;
                    Nodes[i].Y = OutLet.Features.ElementAt(i).Coordinates[0].Y;
                }
                else
                {
                    Nodes[i].X = RandBetween(GridFile.Xllcenter, GridFile.Xllcenter + GridFile.NumColumns*GridFile.CellWidth,i);
                    //Nodes[i].Y = RandBetween(GridFile.Yllcenter,GridFile.Yllcenter + GridFile.NumRows*GridFile.CellHeight,i);
                    Nodes[i].Y = GridFile.Yllcenter + (GridFile.Extent.Height/(double)nNodes)*i;
                }

                dbl[0] = Nodes[i].X;
                dbl[1] = Nodes[i].Y;
                Coordinate cod = new Coordinate(dbl[0], dbl[1]);
                DotSpatial.Topology.Point pt = new DotSpatial.Topology.Point(cod);
                _tree.Insert(dbl, Nodes[i]);
                pt.X = Nodes[i].X;
                pt.Y = Nodes[i].Y;
                DotSpatial.Data.IFeature f = SyntheticPoints.AddFeature(pt);//Add the feature to the FeatureSet
                RI = GridFile.Bounds.ProjToCell(cod);
                try
                {
                    Nodes[i].Z = GridFile.Value[RI.Row, RI.Column];
                }
                catch
                {
                    if (i == 0) Nodes[i].Z = 0.0;
                    else
                    {
                        Nodes[i].Z = Nodes[i - 1].Z;
                    }
                }                          
            }
            MapPointLayer l = new MapPointLayer(SyntheticPoints);
            map1.Layers.Add(l);

        }
        private double RandBetween(double Low, double High, int i)
        {            
            double randomF =(Convert.ToDouble(randonList[i-3].ToString()))/(double)(randonList.Count);                      
            return (High - Low) * randomF + Low;
        }

        private void ProcessSlope(Node CurrentNode)
        {
            double[] d = new double[nNodes - 1];
            int[] keys = new int[nNodes - 1];
            double slope;
            double[] Location = new double[2];
            Location[0] = CurrentNode.X;
            Location[1] = CurrentNode.Y;
            int K = 20; //the number of neighbors;
            object[] Neighbors = new object[K];
            Neighbors = _tree.Nearest(Location, K);
            for (int i = 0; i < K; i++)
            {
                Node Neighbor = new Node();
                Neighbor = (Node)(Neighbors[i]);
                if (Neighbor.Z > CurrentNode.Z)
                {
                    double oldSlope = Neighbor.Slope();
                    slope = Neighbor.GetSlope(CurrentNode);
                    if (slope > oldSlope)
                    {
                        Neighbor.DownNode = CurrentNode;
                    }
                    if (Neighbor.Visited == false) ProcessSlope(Neighbor);
                }
            }
            CurrentNode.Visited = true;
            DotSpatial.Topology.Point HigherPoint = new DotSpatial.Topology.Point();
            HigherPoint.X = CurrentNode.X;
            HigherPoint.Y = CurrentNode.Y;
            HigherPoint.Z = CurrentNode.Z;
            DotSpatial.Data.IFeature f = SlopePoint.AddFeature(HigherPoint);
          }

        private void slopeDelineationToolStripMenuItem_Click(object sender, EventArgs e)//Test the slope based delineation
        {
            ProcessSlope(Nodes[2]);
            DotSpatial.Topology.Point p = new DotSpatial.Topology.Point();
            for (int i = 0; i < Nodes.Length; i++)
            {
                if (Nodes[i].DownNode != null)
                {
                    DotSpatial.Topology.Point[] LinePoints = new DotSpatial.Topology.Point[2];
                    p.X = Nodes[i].X;
                    p.Y = Nodes[i].Y;
                    p.Z = Nodes[i].Z;
                    LinePoints[0] = new DotSpatial.Topology.Point(Nodes[i].X, Nodes[i].Y, Nodes[i].Z);
                    LinePoints[1] = new DotSpatial.Topology.Point(Nodes[i].DownNode.X, Nodes[i].DownNode.Y, Nodes[i].DownNode.Z);
                    DotSpatial.Topology.LineString ls = new LineString(LinePoints);
                    IFeature LineFeature = NetLine.AddFeature(ls);
                }
            }
            MapLineLayer Line = new MapLineLayer(NetLine);
            map1.Layers.Add(Line);
        }
    }
    

    public class Node
    {
        public double X, Y, Z;
        public Node DownNode;
        public bool Visited;
        public double Slope()
        {
            if (DownNode == null)
                return 0;
            else return GetSlope(DownNode);
        }
        public double GetSlope(Node DownNode)
        {
            double dx, dy, dist; ;
            dx = X - DownNode.X;
            dy = Y - DownNode.Y;
            dist = Math.Sqrt(dx * dx + dy * dy);
            return (Z - DownNode.Z) / dist;
        }
    }
}
